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This paper continues development of a full wave propagation theory for low-frequency radio waves 
which is analogous to geometric optics. The effects of earth conductivity, reflection height, and earth 
curvature are described by the path integral, and accurate methods lor computing it, especially suit- 
able for a programmed computer, are given. Knowledge of the path integral and the ionospheric 
reflection coefficient, which are independent, permit calculation of skywave field strengths. Sample 
calculations are shown which confirm that the LF skywave is diffracted deep into the shadow region. 
Possible applications include calculation of LF and VLF field strengths and extraction of ionospheric 
reflection coefficients from field measurements. 

1. Introduction 

The full wave solution for propagation between a spherical earth and a concentric ionosphere 
can be expanded into a series of complex integrals. If the integrals are replaced by their saddle 
point approximations, the series is identified as the ray hop series of geometric optics, so the 
integrals are called wave hops [Berry, 1964]. This expansion was suggested by Rydbeck [1944J 
and Bremmer [1949|. Wait [1961] firmly established the connection with the ray hops, pointed 
out that the saddle point approximation was inadequate near, and beyond, the caustic, and sug- 
gested evaluating the wave hops by numerical integration or by summing a residue series. Wait 
and Conda [1961 J made a few preliminary calculations. Continuing and extending this work, 
Berry [1964] included height gain formulas and calculated the wave hops by numerical integration 
for several representative cases. 

A convenient form of the wave hops series is 

£ r = £o+j7j/i, (1) 

where E is the groundwave and jj is an effective ionospheric reflection coefficient. For a locally 
plane isotropic ionosphere, 

yj = T j , (2) 

where T is the Fresnel reflection coefficient. The effects of ground conductivity, reflection height, 
earth curvature, and distance are accounted for by the path integral, Ij. This paper is concerned 
with calculation of the path integrals. 

Section 2 describes in detail methods for calculating the path integrals for a wide range of 
frequency and distance. Some sample calculations are discussed in section 3, and two possible 
applications are discussed in section 4. 

2. Path Integral Formulas 

A vertical electric dipole source of waves with wave number k = 2iT/\ is on a spherical earth 
of radius a. A time function, exp (iwt), is assumed, and suppressed. The field is to be found on 
the surface at a distance d=a0. The effective reflection height is h. 

For this paper, the starting point is equation (22) of Berry [1964]. Fock's [1946] approximations 
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for the spherical Hankel functions are 



£<«! {x)=-i$)'*W % {t), ( 3a ) 



and 

\ 1/6 /2\ 1/3 



ei w*»^j r ^ t= gj <*-*)■ (3b) 

The Airy functions Wk(x) are discussed by Spies and Wait [1961]. Using (3), and the asymptotic 
form for the Legendre function, P^cos 0), the path integral for the 7th wave hop is 



e -ikd r 

Ij ^ e -isn/4 K —=\ (l+ztfl 2 e- i *W l (t)W 2 (t)(l + R e (t)) 2 Ri- 1 (t)pJ(t)dt, (4) 

Vsin Jr 



(ka\ ^ 3 
where x = vO, y— kh/v, z = 0.5/v 2 , v=\ — ) , 



J) ' and (5) 

T is the contour of integration described below. The ground reflection coefficient is 

KS) - WWWM-q (6) 

where 

q = — v VI — r) 2 jr\ 2 , 

and 

y] 2 = e 2 — ifJLoc 2 (rl(o. (7) 

The relative permittivity of the earth is e 2 and the conductivity is cr; c is the speed of light and /x 
is the permeability of free space. 
Finally, 

P{t) -WMW*t-yY (8) 

and 

K=ll.96hlJ—„v\ (9) 

where 1J is the dipole current moment. In this paper, /(>/= 1. 

Equation (4) differs from that given by Wait [1961] only by the normalization, /£, and the factor 
(l+ztf 12 . Since z < 1, this factor is important only if \t\ is large. A quantitative discussion of 
its effect is given in section 3. 

VL\t\< v\ (1 +ztf> 2 = 1 + 5/2 zt. (10) 

Using the Wronskian [Spies and Wait, 1961], 

W\W 2 -WW x = 2i, (11) 

l + Re[t)= wmw[{!)- q wm (12) 



Then (4) is written 



V^hTeh o +, (t) [10 ' 
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where 

E{t)=W' 2 {t)-qW 2 {t)^ (14) 

F(t)=W 1 (t-y)IW 2 (t-y), and (15) 

C(t)=W[(t)-qW l (t). (16) 

The contour of integration has been discussed by Wait [1961] and Berry [1964]. The contour 
used for the numerical integration in this paper runs along the real axis from oc to the origin, and 
then into the third quadrant on a straight line with slope 1/4. In the third quadrant, the integrand 
gets small exponentially with exp(Im(ta;)), and on the real axis, when t > 1 [Spies and Wait, 1961], 

_J__oc ex p(_|,^, (1?) 

so the integration need be carried out for only a finite portion of the contour in the vicinity of the 
origin. 



2.1. The Saddle Point Approximation 

The development here closely follows Wait [1961]. The difference is a small improvement in 
the accuracy of the approximations. 

U(-t)>\ [Olver, 1954], 

WiAt) = (-0~ ,/4 exp [(-l)*i(f +tt/4)] LK-l)*ifl, (18) 

and 

W' k (t) = (- tyl 4 exp [( - !)*£({ - tt/4)] M [(- !)*#], (19) 



where 



^=|(-«) 3/2 . 

L(Z) = J ^" J . and M(Z) = J FjZ" J , (20) 



where 



Uo=V Q =l, and 

(2/+l)(2/+3)...(6/-l) r;6/±l .., 

" J JK216V • and ^ = - (/j '6M^ >L < 21 ) 

Substituting (18) and (19) into (4), 
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where, letting a = (— 1) 1/2 and s = -a 3 , 



G(t) = (l-ahfl 2 a-\l H-fl) 2 (-^LflY l L\is) 

L{— is) J 



£l(t) = — xa 2 + -j(y + a 2 ) 3/2 — -joe 3 , and 



R 



aM(is) — iqL(is) L(— is) 



aM{- is) + iqL(- is) L(is) 
The saddle point approximation for Ij is (compare Wait [1961]) 



j. a e ~ a ^ 4 K 



Vsin 6 



2n 



iil"(t ) 



1/2 



G(t )e m o\ 



(23) 
(24) 
(25) 



(26) 



where to is the saddle point found from 



Hence, 



lj = 2Ke- inl2 e- ikd 



ao= 4fc £!=( _ fo)1/2> 
4>xj 



1/2 



x sin 6/ \ 2/a , 



L(JSo) 



L(— fco> 



R 



y-i 



L 2 (iso)e 



-ma 



(27) 



(28) 



This reduces to Wait's [1961] result if L, M, and (1 — a 2 z) 5/2 are set equal to one. As a— >°°, 
L and Af— »1; and asa^ 0,(1 — aoz) 5/2 ^T. Thus, the effect of retaining these factors is to increase 
the accuracy of (28) for extreme values of a. Quantitative results are shown in figure 2. 

The reader is referred to Wait [1961] for a complete discussion of (28). He shows that 

x \ 1/2 
1-f— — is the convergence coefficient of geometric optics; H(^ ) is the difference between the 
2ja J 

electric length of the ray path and the distance along the earth; and R is the ground reflection co- 
efficient. Thus, when L = M=1, (28) multiplied by jj is the geometric-optics formula for the 7th 




FIGURE 1. Geometry of first two wave hops. 
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ray hop. This justifies calling yjlj a "wave hop" and suggests using the concepts and terminology 
of geometric optics. Referring to figure 1, <pj is the angle of incidence at the ionosphere of the yth 
hop, and Tj is the angle of incidence at the ground. The ray path length of the /th hop is Dj. The 
region where Tj < 7r/2 is the lit region, Tj = 77/2 is the caustic, and beyond the caustic is the shadow 
reg;ion. 

Equation (27), and the condition ( — t) > > 1, show that (28) is valid only if d/j is sufficiently 
small, i.e., if \tj — tt/2| is sufficiently large. 

2.2. Residue Series for the Path Integrals 
The integrand of (13) has poles of order (j-\- 1) at t s if 

W[(t s )-qW 1 (t,) = 0. ( 29) 

The t s are related to the poles, t s , of the Bremmer [1949] groundwave by [Wait and Conda, 1958| 

t,= ^2T,. (30) 

Wait [1961] shows that these are the only poles, so Ij can be evaluated by summing a series of resi- 
dues—at least in principle. The difficulty is the high order of the poles. Wait [1961] worked out 
the form of the residue for the first wave hop, but did not use it for calculations. Methods for 
calculating the residues of Tj, especially suitable for a digital computer, are given in this section. 
Denote the residue at t s of the integrand of Ij by Res (/, t 8 ). Then 



/, =s 8iriKe i7tl4 



Vsin 6 £ 



£ Res (/, t 8 ). 



(31) 



Rewrite (13) as 






where, using (10), 



Then 



%)a(l+5/2tf)r w 



(33) 





ajo 


Oj\ 


a j2 




bjo 


bp 


bj2 







bjo 


bn 


Res(/, ^) = t Tn? 








bjo 











Ujo 



bjj 

bjj-i 

bjj-2 



bn 



where a Jn (t s ) = Af\t s )/nl, and b Jn (t s ) = B ( j " +j+1 \t,)l(n+j+ 1)1 
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(34) 



(35) 



The superscript (n) denotes the nth derivative. The determinant is easily evaluated numerically 
since it is nearly triangular. 

Using Liebnitz's rule for the nth derivative of a product [Kaplan, 1952], 

~ ml(n-m)l v ; 

the necessary derivatives of [E(t)F(t)] and 

A j (t) = Aj. 1 (t)[E(t)F(t)] (37) 

are found recursively. This repetitive process is particularly appropriate for an electronic com- 
puter. The necessary derivatives of D, E, and F for the first five Ij are 



iy n \t) = (-ix) n - x e- ixt \n-z-ix(l+-zt 



(38) 



E'(t) = tWit)-qW 2 (t), 

E"(t) = tW' 2 {t) + (1 - qt)Wi{t), 
Ef»\t) = (t 2 - q)W 2 (t) + (2 - qt)W' 2 (t), ( 39) 

E (4 \t) = (it - qt 2 )W 2 (t)+{t 2 - 2qW 2 (t), and 
£(s)( ( , = (4 _ 4 gt + t?)JF 2 (t) + (6t - qt 2 )W' 2 {t\ 

F'(t) = 2i[W 2 (t-y)]- 2 , 

F"(t) = -2GF'(t), 

F^={-2(t-y) + 6G 2 }F'(t) (40 ) 

F^(t) ={-2+l6(t- y)G - 24G 3 }F'(t), and 
F^\t)={l6(t-yf + 20G-\20{t-y)G 2 +120G*}F'(t). 



In (40), G= W' 2 (t-y)IW. 2 {t-y). 

The bj„ can be found by repeated use of (36) and (26), but are given explicitly here for the first 
five hops: 

b j0 =(C'(t,)y+\j=l, . . ., 5 b jl =J±l(C'(t s )VC"(t s ),j=l, . . . , 5 

C(t s )b J2 = I C"(t s )b n + ; -ii C*\t s )b }0 , ; = 2,. . . , 5 (41) 

C'(t s )b j3 = J -^ C '\t s )b }2 + ^- C^(t s )b ji + J ~ O*\t,)bjo,j=3, 4, 5 
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C'(t s )b j4 = J -^ C"(t 8 )b J3 + iL C (3) ( ^. 2 + 3Z_^ C ( 4 ) Us)/;ji + ./_M C s) (t8)bjoJ= ^ 5 

6 55 =^ c(t f ) (cw+f (C(t s )f (C"(t s )ro*\ts) 

+ |- (C'(a 4 C (3) (^)C< 4 )(^) + | (CWC'W^ + j^ (C'(* s ))W*s). 

The function C(t) differs from E(t) only by the index on the Airy functions, so using (29) and 
(39) the necessary derivatives of C(t s ) are 

C'(t 8 ) = (ts-q 2 Wi(ts\ 

a*\U) = (%-fU + q)WiU\ 

0*\t 8 ) = {te 8 -2q 2 )Wi{ts\ (42) 

Q 5) (ts) = («J - <7 2 ^ + 2^ + 4)JFi(t f ), and 

C (6) (/,) = (9*J - 6<7 2 /,- + 6q)W x (h). 

Before the convenient form of (34) and the recursive method of (36) and (37) were noticed, 
explicit expressions for Res (/*, t s ),j= 1, . . . , 5, were derived and programmed. To appreciate 
the economy of (34) and (37), note that Res (5, t 8 ) is the sum of 18 distinct terms, and one factor 
of one of these terms, a r >5, is the sum of 73 distinct terms. 



3. Sample Calculations and Discussion 

An efficient program for calculating the path integrals must use all three methods given above. 
The saddle point approximation should be used whenever possible because it is simplest; but it 
is not valid near, or beyond, the caustic. The residue series is accurate and efficient deep in the 
shadow region, but converges very slowly in the lit region. If these two methods do not overlap 
in the caustic region, (13) must be integrated numerically. 

The three methods of calculating the path integrals were compared by computing //, j= 1, 
. . . , 5, as a function of distance, d= 1000, . . . , 8000 km, for 10 kc/s, 100 kc/s, and 1000 kc/s. 
The regions of validity of the methods overlap enough to lend confidence to the computer routines 
and to allow a program to automatically select the appropriate method. Two examples which 
illustrate the comparisons are shown in figure 2. For all calculations in this paper, €2= 15. The 
solid line indicates /,, the other lines show the departure from Ij of the various formulas. The 
simpler saddle point approximation given by Wait [1961] has also been plotted. The caustic is 
indicated by an arrow. Figure 2a shows li (100 kc/s). There is a small difference between the 
value and the simple form of the saddle point approximation near 1000 km. The numerical in- 
tegration also departs from the value here, but this could have been avoided with a higher order 
quadrature. Between 2000 km and 3000 km the error becomes very large in both the saddle point 
approximation and residue series, so numerical integration is necessary. Beyond 6000 km, the 
numerical integration becomes inaccurate. 
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WAVE HOP 4 

a =0.01 mhos/m 
h = 70km 

FREQ = IOkc/S 



SADDLE POINT 

SIMPLER SADDLE POINT 

RESIDUE SERIES 

INTEGRATION 



DISTANCE, KILOMETERS 

FIGURE 2a. Comparison of various methods of calculation, 
h (100 kc/s). 



3000 4000 5000 6000 



DISTANCE, KILOMETERS 



7000 



Figure 2b. 



Comparison of various methods of calculation, 
U (10 kc/s). 



Figure 2b shows I 4 (10 kc/s). For this case, only a small error would result between 2000 km 
and 3000 km from using only the saddle point approximation and the residue series. The simple 
form of the saddle point approximation is not as accurate for 10 kc/s as for 100 kc/s. The error 
is increasing near 1000 km. This is because a is large here, so the approximation, (1 — a 2 z) 5/2 = 1, 
used by Wait [1961], is not valid in this region. However, in most practical applications, the 4th 
hop is insignificant at this distance. 

Graphical presentation of the phase of Ij is impractical because it varies rapidly with distance. 
Instead, a phase lag, relative to the free space ray path, is defined. Let — tt < [A] ^ n denote an 
angle equivalent to A —i.e., [A] differs from A by an integral multiple of 2tt. Then, the path integral 
phase lag, /3, is defined by 

pj = -{ [phase Ij] + [kDj] + tt/2} . (43) 

Since ft(£o) — k(d~Dj) [Wait, 1961], (25) shows that, for perfectly conducting earth (R = 1), [phase Ij] 

IT TT 

= — — — [kd] — [k(Dj — d)]=—[kDj\—— in the saddle point approximation. Thus (3j is generally 

small and varies slowly in the lit region. 

Figures 3, 4, and 5 show the amplitude and phase lag of the first 3 path integrals as a function 
of distance for frequencies from 10 kc/s to 200 kc/s. (Note that the amplitude has been divided 
by the frequency in kilocycles.) For the 70-km reflection height, the first hop caustic is at d = 1880 
km. Figure 3a shows the diffraction of the first hop into the shadow region. As would be expected, 
the lower frequencies diffract further into the shadow region than the higher frequencies. The 
phase lag of I\ is shown in figure 3b. The exponential decrease of the amplitude and the linear 
increase of the phase lag deep in the shadow region confirm the groundwavelike behavior of the 
wave hop beyond the caustic predicted by Wait [1961]. 
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1000 2000 3000 4000 5000 6000 

DISTANCE, KILOMETERS 
Figure 3a. Amplitude of\\ as a function of distance. 



1000 2000 3000 4000 5000 6000 

DISTANCE, KILOMETERS 



7000 



FIGURE 3b. Phase lag of \\ as a function of distance, showing dif- 
fraction into the shadow region. 



Figure 4 shows the amplitude and phase lag of / 2 . The relative minimum in the amplitude 
between 2000 km and 3000 km is caused by the pseudo-Brewster angle in the ground reflection 
coefficient [Bremmer, 1949; Berry, 1964]. The sharpest minimum is for a frequency of 150 kc/s. 
Figure 4b shows the change in the phase lag corresponding to the minimum in the amplitude curve. 
For 150 kc/s there is a rapid phase shift of 180° corresponding to the sharp minimum in the ampli- 
tude curve. The frequency dependence of the pseudo-Brewster angle is clearly shown in figure 
4a. The pseudo-Brewster angle is also shown in the amplitude and phase lag of h in figures 5a 
and 5b. 

The present computer program can calculate h and I 5 but no examples are shown here, since 
no new phenomena occur. Finally, the calculation is not limited to the frequencies shown — a 
test case at 10 Mc/s ran successfully. 

4. Examples of Applications 

Equation (1) shows that the path integrals are the natural complements of ionospheric reflec- 
tion coefficients. Reflection coefficients have been calculated and published for a variety of the- 
oretical models of the ionosphere, l e.g., Johler, Walters, and Harper [1960], Johler and Harper [1962], 
Wait and Walters [1963a, 1963b, 1963c] and Walters and Wait [1963]. Using these, or any other 
reflection coefficients, and graphs of the path integrals, such as figures 3, 4, and 5, theoretical LF 
and VLF field strengths can be calculated quickly and easily by hand. 

For example, Wait and Walters [1963a] show that VLF -reflection coefficients for an expo- 
nential ionosphere are given approximately by |jT| = exp (— A cos <pj) for some constant A. In 
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FIGURE 4a. Amplitude of I 2 as a function of distance, showing the 
pseudo-Brewster angle near 3000 km. 



1000 2000 3000 4000 5000 6000 7000 

DISTANCE, KILOMETERS 

Figure 4b. Phase lag ofh as a function of distance. 



particular, for a frequency of 20 kc/s, and a reasonable choice of parameters, they find A = 3. 
Cos (fj is given exactly by [Bremmer, 1949]: 



COS ifj z 



a 1 1 — cos —\-\-h 



1 - cos — ) - 



-1/2 



(44) 



If jhld«\ 



^ +4 Cir 



(45) 



With this model ionosphere, and the information in figures 3, 4, and 5, the amplitudes of the 
first three wave hops were calculated and are shown in figure 6. The total field could also be 
found by calculating the phases and summing. The sum should include the groundwave at dis- 
tances less than about 1000 km. 

The path integrals can also be used to remove the effects of the earth's conductivity, diffrac- 
tion loss, etc., from experimental data. In figure 7, monthly averages of 100-kc/s first-hop sky- 
wave amplitudes are shown as a function of GMT for a 2510-km path from Attu to Sitkinak in the 
Aleutians [Doherty, 1964]. These are Loran-C pulse measurements, so the groundwave could 
be measured before the first-hop skywave arrived. The ordinate scale is A = 20 logio (EiIEq)+ 12. 
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Figure 5a. Amplitude of hi as a function of distance. 



2000 3000 4000 5000 6000 7000 
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FIGURE 5b. Phase lag of I 3 as a function of distance. 



Since E\ =77i, the ionospheric reflection coefficient is given by 



log™ T= 



A-Yl 
20 



Ogio 



h 

Eo 



(46) 



The ratio of the first path integral to the groundwave was calculated for reflection heights of 
65 km (daytime reflection height) and 85 km (nighttime reflection height), for a sea path. These 
ratios are 217.4 and 295.3, respectively. With these ratios, and (46), the measurements in figure 7 
can be converted directly to reflection coefficient magnitudes, without making assumptions about 
the antenna pattern or transmitter power. Note, however, that the reflection coefficient is depend- 
ent on the assumed reflection height and the assumed ground constants. For example, in Sep- 
tember 1962, the average grazing reflection coefficient during the day for this path ranged from 
0.013 to 0.018, and was about 0.07 for several hours during the night. 

5. Concluding Remarks 

The major barrier to practical use of the wave hop theory has been the difficulty of calculating 
the path integrals near the caustic and in the shadow region. The necessary formulas for such 
calculations are given in this paper and sample calculations are shown illustrating the general 
characteristics of the path integrals. The path integrals can be used in conjunction with pub- 
lished values of reflection coefficients to calculate theoretical LF, VLF fields by hand, or they 
can be used to extract reflection coefficients from experimental data. Therefore, the path integrals 
are being calculated as a function of distance for seven frequencies from 10 kc/s to 200 kc/s, for 
reflection heights from 60 km to 100 km, for sea water paths, and for land paths with conductivity 
0.01 and 0.001 mhos/m, and will be published soon [Berry and Chrisman, 1965J. 
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TYPICAL AMPLITUDE CURVES FOR SITKINAK VERSUS 
HOUR OF DAY AND VERSUS ZENITH ANGLE OF THE SUN 



DISTANCE, KILOMETERS 



Figure 6. Sample calculation of amplitude of first three hops, 
using figures 3, 4, and 5, f = 20 kc/s, |T| =exp {—3 cos <p). 




FIGURE 7. Measured first-hop skywaves relative to groundwave, 
{=100 kc/s [after Doherty, 1965]. 
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